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ABSTRACT 

Context. We have investigated the structure and kinematics of triaxial final models formed due to radial-orbit instability from a set of 
equilibrium anisotropic spherical systems of the Osipkov-Merritt type. 

Aims. We show that the instability is a natural way to build nonspherical systems and thus can be considered as a new technique for 
constructing equilibrium ,/V-body models. 

Methods. The dynamical evolution of the models was followed numerically by means of the Dehnen's public code gyrfacON. 
Results. We found that the shape of the density profiles for various initial spherical models didn't change despite of the systems 
were drastically rearranged into triaxial configurations and had got the new scale lengths. The cusp was found to be an invariant 
feature of such a rearranging. Starting from a certain unstable spherical y-model we obtain a triaxial system with the same cusp. The 
end-products are anisotropic at large radii and have anisotropy profiles of the Osipkov-Merrit type. The size of the isotropic core once 
rescaled corresponds to the value of marginally stable y-progenitors and is almost independent of the unstable starting point. 
Conclusions. We conclude that the end-products reach a new steady state and have quite predictable properties. They can be used 
as equilibrium models to describe the elliptical galaxies and bulges. They also can be incorporated in ,/V-body simulations deal with 
multicomponent systems. 

Key words, stellar dynamics — instabilities — methods: A'-body simulations — galaxies: structure — galaxies: kinematics and 
dynamics 



1. Introduction 

Triaxiality is very widespread among stellar systems. For ex- 
ample, there are a few perfectly spherical elliptical galax- 
ies. The analysis of their morphology and kinematics have 
shown that elliptic als are at least moderately triaxial systems 
dFranx et al.l |1991[). Dark matter haloes formed in isolation 
(e.g.. lAguilar & Merritt &1 Il990t iDubinski & Carlberd [l99lt 
lKat3l799lT or in subregion s of larger cosmolo gical simula- 
tions (e.g., lFrenketalJll988t IWarren et allll992l) through dis- 
sipationless collapse of density peaks also demonstrate triax- 
ial shapes. The same picture arises in cosmolog ical simulations 
with halo assembling via merging (see, e.g.Jjing & Suto| 2002; 
Bailin & Steinmetd2005tlShaw et all2006tlAllgood et al.ll2007t 
Bett et alj|2007l) . 



There is one more common feature of real and cosmologi- 
cally motivated systems. It is a lack of constant density cores. 
As a rule such systems have cuspy density profiles. The stellar 
surface brigthness of early-type galaxies very often continues 
to rise towards small di stances from the centre (Ferraras e et al.l 
1994t lLauereTaT] [T995). But even a core-like surface brightness 
distribution can be transforme d due to projection into a power- 
law cusp in its spatial density (lLauer et al ] |1995l) . As is known 
for dark matter haloes those have universal cuspy density dis- 
tributions — the so- called NFW profiles with p oc r~ y near the 
centre, where y — 1 dNavarro et al ] |1996Lll997h . 

It is very important to have different reliable prescriptions for 
constructing equilibrium models of such systems and to know 



their structural and dynamical properties. There are at least two 
problems which require such models. One of them deals with 
photometric and kinematic data fitting. The other includes many 
question connected with dynamical evolution of triaxial stellar 
systems in a controlled manner and related topics (e.g., the ef- 
fects of triaxial halo shape on the disk stability and bar evolution, 
the survival of triaxiality during the disk formation and so on). 
The most practical models to solve the later tasks are A^-body 
realizations of stellar systems. There are various and well devel- 
opped techniques for constructing equilibrium models, mainly 
invented for spherical ones with known DF (for the DF con- 
struction of models with given density and v elocity anisotropy 
profiles see IOsipkovlll979l: lMerritdll985a|lbl; ICuddefordlll99li: 
iBaes & Deionghdl2002t iBaes & Van Hesdl2007l) . or for the su- 
perpo sition of spherical (or spheroid) and disk components 
(e.g. iHernquistl Il993b iBoily et al.l |200U iKuiiken & Dubinsk " 



ski 
)7; 



1995; Widrow & Dubinskjl 119951: McMillan & Dehnenl 1200 
Rodiono vet al.ll2008l) . 

To build triaxial models with given structural and kinematic 
parameters is a far more complicated problem. Jeans theorem 
and a knowledge of the DF is useful only in a few special an- 
alytical cases. The most popular numerical tech nique for triax- 
ial sy stem modelling is Schwarzschild's method dSchwarzschil 
1 1979 ). Tested initia lly on core-like systems (Sch warzschiL 
Il979[ IStatlerl 1 19871) it was further applied to obtain self- 
consistent gal axy models with cusps ( iMerritt & Fridm an 1996; 
lMerritt|[l997l) and recently has been generalized to construct 
models of triaxial elliptical galaxies embedded in triaxial dark 
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matter haloes dCapuzzo-Dolcetta et al.l 120071) . These models 
based on a library of orbits in a given potential are well suited 
to investigate the variety of orbit types presented in a model but 
have almost never been used to produce initial conditions for 
Af-body simulations so far. 

iHollev -Bockelmann et al] d200ll) proposed a technique 
called 'adiabatic squeezing' to produce cuspy triaxial A^-body 
models. They st arted with a spherical isotropic Hernquist model 
dHernquisll 19901) and eventually were applying a drag to the par- 
ticle velocities along each principal axis. During the squeezing, 
the system retained its initial cuspy density profile with y — 1 
and displayed a slight radial anisotropy in the velocity ellipsoid. 
Rodi onov et al.l (120081) described a very general prescription to 
construct a triaxial model by means of an iterative metho d (in 
original version proposed by iRodionov & Sotnikovall200"6h and 
gave an example of such a model with a given density distribu- 
tion and several constrains on kinematics. 

All these techniques give us equilibrium or near equilibrium 
triaxial models but tell nothing about the nature of triaxiality. 
Triaxiality can arise in a number of ways. One natural way seen 
in hierarchical cosmological simu lations is the merger of a num- 
ber of objects. Moo re et al.l (12004) explored the generation of tri- 
axial structures formed via merging isotropic equilibrium spher- 
ical haloes of NFW type with varying amounts of angular mo- 
mentum. They investigated the relationship between the value of 
angular momentum and the final shape of configurations (prolate 
or oblate). Despite of the final shape their resulting haloes had 
the same density profile as the progenitor haloes, independent of 
their angular momentum. 

Another common way of triaxial system formation is in- 
stability. Starting from equilibrium but unstable conditions a 
stellar system can achive a new equilibrium that can be sta- 
ble. The leading instability pushing the system towards a 
new equilibrium and rearranging its initial spherical density 
distribution into a triaxial one is the radial-orbit instability 
(Polv achenko & Shukh manlll981l: for a physical description of 
this instability see Fridman & Polvachenko 1984). It manifests 
itself in systems with a high degree of velocity anisotropy 
dMerritt & Aguilarll 19851: iMeza & Zamoranolll997t iBuvle et al.l 
12001) . This instability is tho ught to be the main driver for the for- 
mation of elliptical galaxies dvan Albadall982l) . Sufficiently cold 
initial conditions force system stars to fall almost radially toward 
the centre. And if the system has too many radial orbits such 
a collapse will result in a t riaxial configuration. Motivated by 
these ideas iTrenti & Bertinl (|2006) fulfilled some numerical ex- 
periments to clarify the following questions: what kind of equi- 
libria can be achieved in the process of the collisionless collapse 
starting from various conditions, what are the structural proper- 
ties of the end-products and are they similar to the real galaxies. 

In this paper we address to initially equilibrium spheri- 
cal models with cusps and different degrees of the velocity 
anisotropy and analyse shapes, structure and kinematics of final 
triaxial configurations produced via radial-orbit instability. 

The outline of this paper is as follows. In Section 2, we de- 
scribe the technique for constructing initial equilibrium but non- 
stable cuspy A^-body models with various degrees of velocity 
anisotropy and the A^-body code used to study their evolution. In 
Section 3, we investigate the dynamical evolution of our models 
towards a new equilibrium and analyse their structure and kine- 
matics. In Section 4, we discuss the properties of the final models 
in context of application to real stellar systems. In Section 5, we 
summarize our main results. 



2. Modelling technique 

2.1. Family of initial models 

We start with cusp y spherical systems from the one-p arameter 
family y-models (Dehnenl ll993l; iTremaine et alj 1 1994b . These 
models follow a density profile given by 



P(r) = 



(3-r)fl 



M 



An r y(r + df-y ' 



(1) 



where M — the total model mass and a — the typical scale 
length. The density of the models behaves asp oc r~ 4 for large 
radii and as p cc r~ 7 near the centre, thus y characterizes how a 
central cusp is strong. 

We have limited our study to models with y — 1 ('weak- 
cusp') and y = 3/2 ('mode rate-cusp'). The former model is the 
so-called Hernquist sph ere (|Hernciuist| [i990). Surface density of 
both models is shown (Dehnen 1993) to agree q uite well with 
the de Vaucoulers profile (de Vaucouleurs 1948) that describes 
the surface brightness of many ellipticals and bulges of spirals. 

We consi dered a set of anisotropic mod els of the Osipkov- 
Merritt type dOsipkovlll979l: lMerrittl[T985allbl) . In this case the 
anisotropy parameter f5 — 1 - crf/2cr%, where cr\ and erf are the 
tangential and radial compone nts of the velocity dispersion, is a 
function of radius in the form dOsipkovll 1 9791; iMerritl 1985allbh 



OM 



r z + r: 



2 ' 



(2) 



The anisotropy radius r a cut off the isotropic core from the 
radial-orbit dominated region outside r a . 

Models with a given density distribution have two critical 
values for r a . One value seperates physical models with a strictly 
positive DF fro m unphysical ones. Anothe r value marks the sta- 
bility b oundary. iMerritt & A guilar ( 1 9851) : IMeza & Z amorano 
dl997l) : lBuvle et al.l d2007l) set numerically the stability thresh- 
old for various models including y-models of Osipkov-Merritt 
type. For y = 1 the critical value of r a is w 1.1, for y = 
3/2 the stability thresho ld was found to be at r a « 0.8 
(M eza & Zamorano! 1 19971) . We consider long time dynamical 
evolution of physical models lying below the stability threshold. 

In all numerical simulations we used the following dimen- 
sionless units: G = \, M = \, a = 1. These units can be con- 
verted in physical ones through typical observational parame- 
ters of elliptical galaxies — mass and effective radius that char- 
acterizes the d e Vaucoulers pr ofile. For y — 1 effective radius 
r e = 1.815a dHernquistl 119901 for y = 3/2 it can be derived 
from the cumulative surface density dDehnenl 1993b r e — 1 .276a. 
If we adopt M — 2 • 10 1 l M®, r e — 3 kpc as typical parameters for 
bright ellipticals it will give us the scale unit r„ = a = 1 .65 kpc, 
the time unit t„ = 2.24 Myr, the velocity unit v„ = 721.3 km/s 
fory = l,andr„ = a = 2.35kpc,f„ = 3.80Myr, v„ = 604.7km/s 
for y = 3/2. 

IMeza & Zamorano! (fT997h : lBuvle et alJd2007l) expressed the 
evolution of their models in dimensionless units of half-mass 
dynamical time T^. It is ev aluated at the half-mass radius r\ji 
dBinnev & Tremainelll987l) . For y = 1 it is equal 7h ~ 8.3 
and T h ss 4.9 for y = 3/2. lDehnenl d 19931) noticed that the ra- 
tio r e /ri/2 ~ 0.75 is practically the same for all y-models. This 
means that if we fix the physical parameters (mass and effective 
radius) we will obtain the same physical dynamical time scale 
for all models. For y = 1 with adopted parameters it will be 
T h = 18.7 Myr and T h = 18.8 Myr for y = 3/2. 

We followed the dynamical e v olution of the mod els further 
than in lMeza & Zamoranol d 1997b : lBuvle et all (120071) . for some 
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models - up to t * 1000 - 1400 (i.e. * (100 - 200)r h ). It con- 
cerns the models just below the stability threshold. They were 
approaching the new steady state very slowly exhibiting a down- 
fall to a new constant level in the behaviour of the axis ratio only 
after 70 - 120 dynamical times. 

2.2. Initial data sets 

An iV-body realization of an initialy spherically symmetric sys- 
tem was created in a conventianal way by using the DF. For 
Osipkov-Merritt models the DF has a form 



temporal behaviour of the axis ratios b/a and c/a as indi- 
cators of new steady state settling. Till now such investiga- 
tions were concentrated on finding numeri cally the stability 



tions were concentrated on finding numerical ly me sta bility 
threshold of spherical anisotropic models (|Merritt & Ag uilar 
Deionghe & Merrittl Il988t IMeza & Zamoranol 119971 



1985 



Mezall2002t Buvl e et al. " 120071) . Triaxial end-products formed 
due to the radial-orbit instability from initially equilibrium sys- 
tems were commonly considered as a 'garbage'. Below we stud- 
ied the structure and kinematics of such systems settled in a new 
equilibrium in details and determined the relationship between 
properties of initial and final models. 



f(S,L)=f(Q), 



(3) 



where Q — S - L /2r\, & — *F - jv , *F is the negative of 
the gravitational potential, and L = r\\ is the angular momen- 
tum with v t the tangential velocity components. The DF in the 
form (O f or a given density profile can be found through the 
technique (Osirjkov ll979t lMerritl|[T985 a.b) which is an exten- 
sion ofEddington'sJmS) inversion technique to obtain the DF 
for an isotropic spherical system. 

The initial conditions in the A^-body problem suggest speci- 
fying the mass, position in space, and three velocity components 



0v 



Vg) for each particle. 



The mass of all particles was assumed to be the same. As the 
DF is specified the particle positions and velocities are initial- 
ized by sampling the DF. The particle coordinates are naturally 
determined in accordance with the density profile (fTJ. The usual 
way to do this is to invert the cumulative mass profile for the 
radius. 

For velocity components it is enough to sample v r and 

Vt = ytfp + v e an d to cnoose the random direction for v t in 
the plane which is perpendicular to the radius-vector. The algo- 
rithm for velocity sampling in anisotropic systems of more gen- 
eral type than the Osipkov-M erritt models is breifly described in 
McMillan & Dehnen (2007) and implemented in the public code 
mkhalo which we have taken from the Walter Dehnen' s part of 
the p ublic NEMO package (http://astro.udm.edu/nemo; iTeubenl 

2.3. Computational method 

We investigated the dynamical evolution of our mode ls by using 
the N-body code gyrfalcON dDehnenll2000l 12002b that com- 
bines a hierarchical tree method (Barnes & Hutll 19861) and a fast 
multipole method dGreengard & Rokhlinlll987l) . This code has 
complexity O(N) and is faster than a standard tree code by at 
least a factor 10. Its impleme ntation was taken from the public 
NEMO package (T eubenll 19951) . 

For almost all models we chose the number of particles N = 
300000. The tolerance parameter 6 that is responsible for the 
accuracy of calculating the gravitational force was fixed at 0.6 
in all our simulations. Calculations were performed in a single- 
time-step mode. In this case momentum is exactly conserved. 
The time step and softening length were taken dt - 1/2 6 and 
e = 0.04. The choice of these parameters was in agree ment with 
the recommendations of Rodi onov & Sotnikoval y005). 



3. Properties of final equilibrium models 

We performed the simulations on a long time scale (up to 
about 200 dynamical time for some models), and used the 



3.1. Model shapes 

To determine the shape of the particle distribution s we used 
a tech nique desci bed bvlDubinski & Carlberd dl 99 lb andlKatzl 
d!99ll) (see also IMeza & Zamoranol 1 1997b iBuvle et alj l2007lh 
First of all we transformed the input A^-body snapshot at spec- 
ified time moments to a coordinate system which diagonalizes 
the moment of inertia tensor /y = 2 XiXj/r 2 of a specified sub- 
set of particles within a certain sphere. For the sphere radius we 
chose r = 5. The centroid of the subset was translated to the 
origin. After diagonalization the principal components of the in- 
ertia tensor I xx > I yy > I zz were used to calculate the preliminary 
values of the axis ratios 



a Uxx a \I XX 



1/2 



(4) 



After that we determined the ellipsoidal radius q as q = (x 2 + 
y 2 /(b/a) 2 + z 2 /(c/a) 2 ) 1 / 2 and calculated the eigenvalues and 
eigenvectors of inertia tensor for a subset of particles with q < 5. 
Tensor components were found as 7, ; = T,(XiXj/q 2 ). The new 
axis ratios were used as the conditions for the next iteration. We 
repeated this procedure until the error of both axial ratio deter- 
mination became less than 10~ 4 . 

Figures[TJand|2]demonstrate the temporal behavior of the ra- 
tios between the intermediate and major axes and minor and ma- 
jor axes for two y-models with different values of the anisotropy 
radius r a . 

Reaching the constant value of the axes ratio was consid- 
ered as a new equilibrium. Some models with the anisotropy 
radius just below the stability threshold were approaching this 
state on a time scale longer than 7h = 50. As a result we ob- 
tained more extended end-products. F or example, for y = 1 the 
stability threshold was set as r a = 1.1 ([Meza & Zamorano 1997; 
IBuvle et aUl2007l) . Starting from r a = 1.05 we have got after 
150 dynamical times (f « 1200) a perfectly prolate model with 
b/a, c/a 0.77. The model with the anisotropy radius far from 
the stability threshold (r a ) seems to be in nonequilibrium even 
after 100 dynamical times. A slow axes evolution is not notice- 
ble on a short time scale (7h = 50) but is clearly seen during a 
longer period. 

In the Table [TJ we summarize the properties of all sets of 
simulations. The values of b/a and c/a are averaged over last 
100 units of the time. They are very close to each other. It means 
that all the end-products of the initially unstable y-models have 
an extremly prolate shape with a slight triaxiality for the model 
with y = 1 and r a = 0.5 and for the nonequilibrium model with 
y = 1 and r a = 0.3. 
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Fig. 1. Longtime evolution of the intermediate b and minor c to major a axis ratios for the y — 1 models with different values of 
initial anisotropy radius r a . The axes lengths are iteratively calculated from the ellipsoidal density distribution using the moment of 
inertia tensor. 
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Fig. 2. Longtime evolution of the b/a and c/a ratios for the y = 
3/2 models with two values of initial anisotropy radius r a . The 
axes length are iteratively calculated from the ellipsoidal density 
distribution using the moment of inertia tensor. 



Table 1. Properties of the final models 



y 


r a 


'end 
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a 


b/c 


c/a 




l 


1.05 


1400 


1.46 


2.11 


1.11 


0.78 


0.77 


1.31 


l 


0.9 


1200 


1.59 


2.07 


1.24 


0.65 


0.65 


1.29 


l 


0.7 


1000 


1.90 


1.94 


1.53 


0.51 


0.51 


1.24 


l 


0.6 


1000 


2.00 


1.79 


1.83 


0.47 


0.46 


1.10 


l 


0.5 


400 


2.08 


1.79 


1.93 


0.45 


0.43 


1.07 


l 


0.3 


800 


1.91 


1.64 


2.03 


0.48 


0.39 


0.94 


3/2 


0.7 


1000 


1.21 


2.02 


1.06 


0.75 


0.75 


0.92 


3/2 


0.5 


1000 


0.97 


2.09 


1.44 


0.53 


0.53 


0.84 



3.2. Density profiles 

The most remarkable feature of the end-products is their density 
profiles. During reaching a new equilibrium the y-models retain 
their initial density profiles. Once the new scale length a is found 



all density profiles of the final models can be described by a law 
likes (Q]i but modified in the case of triaxiality as 

n _ (3-y)a 3 l 

P q 4nbc qy{q + a) 4 ~v ' 

where q is ellipsoidal radius. 

Figure [3] shows the density of the end-products as a func- 
tion of an ellipsoidal radius q. We demonstrate also the fitting 
of some final density profiles by formula © with b/a and c/a 
taken at the end of simulations. The profiles can be fitted with an 
accuracy of better than 1%. The models exhibit an unchanging 
density profiles over ~ 3 - 4 orders of magnitude in radius af- 
ter settling to the new equilibrium, and they do not evolve away 
from the original cusps. The density slope seems to be an invari- 
ant feature of these models despite of their drastical rearranging 
into the triaxial configurations. Starting from a spherical model 
with a certain y cusp we obtain a triaxial model with the cusp of 
the same shape. 

3.3. Kinematics and velocity anisotropy 

The velocity dispersions (cr,., cr^, erg) as a function of radius 
in Osipkov-Merritt models can be easy compu ted as the cor re- 
sponding moments of the DF in the form f(Q) (lMerrittll985allbh . 
The y-models of the Osipkov-Merrit type have a distinctive 
hump in their velocity dispersion profiles provided y < 2. 
The same feature manifests its elf in the radial pro files of the 
projected velocity dispersion (Caroll o et al.l 1 19951) . The end- 
products exhibit the same hump in their velocity dispersion pro- 
files like their spherical progenitors (Fig. 0). Despite of their 
nonspherical shape all models demonstrate the equivalency be- 
tween v v and Vfl velocity components as the spherical models 
with DF in the form /(£, L). 

All final models display a radial anisotropy in the velocity 
ellipsoid. We choose to represent the anisotropy profiles with 
more general law than (fj) 

with r a f and 6 being free parameters. The results of fitting are 
summarized in the Table Q] All models obtained from progeni- 
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Fig. 3. Density of the models plotted as a function of ellipsoidal radius q. The curves show the density profiles for the final triaxial 
models formed via the radial-orbit instability from the initial y-models with different values of anisotropy radius r a . The thin solid 
lines show the fitting of profiles by formula (O. Curves are shifted from one another by dex=0.5 (from top to botom) for convenience. 



tors laying just below the stability threshold have the anisotropy 
profile with 6 ~ 2, e.g. they are of Osipkov-Merritt type. Fig [5] 
shows some examples. The profiles were fitted over the range 
q < 5. The deviation from (|6]l are observed only at large dis- 
tances q > 5 from the center and are due to a slight predom- 
inance of tangential orbits there as compared with Osipkov- 
Merritt models. 



4. Discussion 

(1) Triaxial models formed due to instabilities from spherical 
anisotropic equilibrium systems were commonly thought as a 
'garbage' . Their structure and kinematics weren't investigated in 
details. We showed that such models can be considered as equi- 
librium ones and used for studying their dynamical evolution 
and related topics concerned with dynamics of multicomponent 
systems. 

(2) All our models have a near prolate shape without any 
strong triaxiality. Such a shape is a common feature of cosmo- 
logically motivated dark halos. There have been many theoret- 
ical papers published over the years which examined the sub- 
ject of halo shapes. Most authors found that halos tend to be 
prolate. For example, using six high resolution dissipationless 
simulations with a varying box size in a flat LCDM universe 
Allg ood et al.l d2007l) found that most halos in their simulations 



were prolate in shape with very few oblate halos. But our pro- 
late models, that are based on initial spherically symmetric y- 
models, are rather suit to describing the density distribution of 
ellipticals and bulges. 

Several authors discusse d the problem of the intrinsic shape 
of elliptical galaxies (e.g., | Franx et al.l 1199 it IVincent & Rvdenl 
l2005t iPadilla & Straussll2008l) . Many of them agree that the ob- 
served distribution of ellipticities can be reproduced by suggest- 
ing the triaxiality of these objects. It is concluded that any ac- 
ceptable distribution is dominated by nearly -oblate spheroidal 
rather than nearly-prolate spheroidal systems. Padilla & Straussl 
(2008) have recently reexamined the underlying shapes of el- 
liptical galaxies in the SDSS Data Release 6 from the observed 
distribution of projected galaxy shapes, taking into account the 
effects of dust extinction and reddening. The elliptical galaxy 
data were found to be consistent with oblate spheroids, with a 
correlation between luminosity and ellipticity: more luminous 
ellipticals tend to be rounder, although ellipticals are oblate at 
all luminosities. Thi s resul t is in agreement with an analysis 
bv lVincent & Rvdenl d2005l) based on SDSS Data Release 3 al- 
though later authors found that the fainter de Vaucoulers galax- 
ies are best fitted with prolate spheroids with mean axis ratio 
<(A,c)/a>!»0.51. 

iMendez-Ab reu et all (120081) performed two-dimensional 
photometric decomposition of the galaxy surface brightness dis- 
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Fig. 5. The velocity anisotropy parameter as a function of ellipsoidal radius q. Thin curves — the initial profiles, thick curves — the 
profiles for end-products, thin dotted lines show profiles obtained by fitting the data by formula @. Left — 7=1; solid lines — 
r a = 0.9, dashed lines r a = 1.05. Right — y — 3/2, solid lines — r a = 0.5, dashed lines r a = 0.7. 



tribution to derive the structural parameters of disks and bulges 
for 148 unbarred SO-Sb galaxies. They obtained the probabil- 
ity distribution function of the intrinsic equatorial ellipticity of 
bulges and found that about 80% of bulges in their sample are 
not oblate but triaxial spheroids. This is consistent with several 
previous findings. 

Thus, in the context of observational data our models can be 
related only with the lowest luminosity ellipticals. 

(3) All our final models are strongly elongated with 
{(b,c)/a) ~ 0.5 - 0.8. It's true even for the models start- 
ing from the conditions just below the level of stability. For 
r a = 1.05 we folowed the dynamical evol ution of the system 
further (T en d ~ 1707a) t han it was done by iMeza & Zamoranol 
d 1997b and iBuvle et al.1 (120071) and obtained a very elongated 
end-product. 

(4) One of the main points of our investigation is that the 
shape of the density profile didn't change. The end-product 
structure can be predicted in a qu ite definite w a y and depends 
only on initial model properties. iMoore et al.l (120041) noticed 
that their triaxial remants formed via merging keep cen tral den- 
sity slopes of progenitors. iBovlan-K olchin & Ma (2004) claimed 



that if the cuspy halo had started out with an NFW form with 
7=1, the merged halo also had a form close to NFW one, 
with a slightly shallow in n er cus p of p oc r _0i instead of r . 
iHollev-Bo ckelmann et al. (2001), who had used a squeezing 
technique for constucting triaxial systems, also concluded that 
their models retained the given y character after squeezing had 
terminated. It seems the cusp shape is an invariant feature of any 
matter rearranging as though the particles retain the memory of 
the progenitor density profile. 

(5) The dark halos obtained in cosmological simula- 
tions have a small degree of anisotropy far from the center 
(ICole & Lacevl Il996t ICoim et al.1 120001 iFukushige & Makinol 
120011: iDiemand et al.ll2004b . We know a little about 3D veloc- 
ity distributions in ellipticals. But modern methods for mod- 
elling of observational data allow to reconstruct such distribu- 
tions. Combining surface bri ghtness photometry and long-slit 
absorbtion line kinematic data. lDe Lorenzi et all (2008) have re- 
cently presented a dynamical study of NGC 4697, intermediate- 
luminosity, E4 galaxy. Their best fitting models were only 
slightly radially anisotropic, with the anisotropy parameter in- 
creasing to j3 0.5 at large radii. 



T. Piffl and N.Ya. Sotnikova: The properties of cuspy triaxial systems formed due to the radial-orbit instability 



7 



The kinematic properties of the prolate models under study 
aren't similar to the properties of the cosmological dark ha- 
los and ellipticals. The possible explanation is that we started 
out from very anisotropic models of the OM type. The velocity 
anisotropy profile on the periphery was produced by purely ra- 
dial orbits (J3 tends to 1). The end-products are also anisotropic 
at large radii but [} is less than 1 (see Fig [5]). However, as one 
can see from the Table [1] the size of the isotropic core increased 
for all models and once rescaled it corresponds to the value for 
marginally stable y-progenitors and almost independant on a 
starting unstable point. As opposed to this picture, in a merging 
scenario the particle velocities, being initially isotropic, remain 
mostly isotropic near the centre of the merger remnant, but be- 
come mildly anisotropic (fi < 0.4) with inc reasing radius due to 
radial infalls (iBovlan-Kolchin & Mdl2004l) . 

Velocity dispersion profiles of our models demonstrate one 
remarkabe feature — the equivalency between and vg velocity 
components. This fact in combination with a model shape may 
help to construct the DF of such prolate systems in the form 
/(£, L). 

(6) We followed the evolution of our models further than 
in pre vious investigations dMeza & Zamoranoll997tlBuvle et all 
120071) . These authors studied the stability of their models by us- 
ing an jV-body code based o n the 'self-consistent field' method 
dHern quist & Ostrikeri ri992l) . in which the density and the grav- 
itational potential are expanded in a biorthogonal set of some 
basis functions. The using of spherical harmon ics provided a 
reflection symmetry . We used a tree-like code (|Teubenlll995t 
lDehnenll200d 12002) and reached a definitive equilibrium only 
by a special symmetrization procedure that preserves a reflec- 
tion symmetry of a model during a simulation (the manipulator 
symmetrize in Dehnen's gyrfalcON). One have to keep this 
fact in mind when modelling such systems. 
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Fig. 6. Longtime evolution of the c/a ratio for the y — 1 models 
with r a = 0.7. Thin curves — calculations without symmetriza- 
tion, thick curves — calculations with symmetrize manipula- 
tor. The axes length are iteratively calculated from the ellipsoidal 
density distribution using the moment of inertia tensor. 

(7) lMerritt & Fridmanl (fl996l) : lMerrittl d 1997b considered the 
question about the degree to which triaxiality can be supported 
in galaxies with high central densities. They concluded that tri- 
axiality is inconsistent with a strong cusp due to rapid chaotic 
mixing and a real galaxy might choose axisymmetry with a ten- 
dency to be nearly prolate spheroidal, nearly prolate, or nearly 



spherical. All our models formed via the radial orbit insta- 
bility are axissymmetric with a prolate shape. This form is 
reached on a time scale of about ~ 50 7^. After that it can be 
preserved only by a 'hand-made' maintenance of a reflection 
symmetry. Without this maintenance some of our systems un- 
dergo the instability and start evolving towards nearly oblate 
spheroids (see Figure |6). It is known that instability criteria are 
the strongest constraints on structural and kinematic parameters 
of stellar systems and instability is a natural way for equilibrium 
stellar system formation. That is why the models constructed 
here need further investigations. Futher studies will analyse the 
long time orbital evolution of particles to study orbit families 
consistuting these models and the degree of chaos. 

5. Summary 

Our results can be summarized as follows: 

1 . The radial orbit instability is a natural way to build nonspher- 
ical systems and can be considered as a new technique for 
constructing equilibrium A^-body models. 

2. This models have quite predictable properties because the 
shape of the density profiles doesn't change in the process of 
instability. For example, the initial cusp survives. 

3. New steady state models have a definitely prolate shape. 

4. The degree of anisotropy is diminishing and end-products 
are more isotropic. 
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